clear;clc;
addpath('fun')

rng('default')

% parameter values (annual)
s     = 2;          % IES (standard, Zhou)
beta  = 0.95;      % discount factor (BFHV)  
theta = 0.2;        % down paynebt ratio (standard)     
ra    = 0.01;       % interest rate on liquid assets (data Zhou)
fb    = 0.03;     % refi closing cost

rb_bf = 0.04;   % ss mortgage rate (data, Zhou)
rb_af = 0.03; 

mar   = 0.001;       % prevent c=0
infi  = 1000;        % prevent c=0, or interp out of bound
diffini = 100;
difftol =.0001;

alpha=3;

  % pack into a parameter vector
par{1}=s; par{2}=beta; par{3}=theta; par{4}=ra; par{5}=fb;
par{6}=mar; par{7}=infi;par{8}=diffini;par{9}=difftol;par{12}=rb_af;par{13}=alpha;


%set up grids
amax = 1;
amin = -0.01;

na  = 16; 
na1 = na; 
na1_opt = 2*na1;
nb1 = 6;

% endogenous states (h,b,a)
amat       = getnodes(na,amin,amax);
globvar{1} = amat;

% log-income process:z=log(y);z(t)=(1-rhoz)muz+rhoz*z(t-1)+ez;ez~N(0,sez^2)
nz       = 7; 
mz       = 2;  
rhoz     = 0.9;
sez      = 0.1;      
sz       = sez/sqrt(1-rhoz^2);
muz      = -sz^2/2;   % for exp(z)=1, log normal mean

% Tauchen(1986)- discretize AR(1)
[zmat,ztran] = tauchen(muz,rhoz,sez,nz,mz);

globvar{2}   = zmat;
globvar{3}   = ztran;

% log-hp process:p=log(P);p(t)=(1-rhop)mup+rhop*p(t-1)+ep;ep~N(0,sep^2)
pyratio  = 4; 
np       = 7; 
mp       = 2;  
rhop     = 0.95;
sep      = 0.05;      
sp       = sep/sqrt(1-rhop^2);
mup      = log(pyratio)-sp^2/2;   % for exp(z)=1, log normal mean

% Tauchen(1986)- discretize AR(1)
[pmat,ptran] = tauchen(mup,rhop,sep,np,mp);
globvar{4}   = pmat;
globvar{5}   = ptran;

bmat  = (1-theta)*exp(pmat);

% Baseline utility cost
fH    = -1.25; 
par{10} = fH;

pf    = 0.75; 
fmat  = [fH 0]'; fprob1= [pf 1-pf]; nf = length(fmat);
fprob = repmat(fprob1,nf,1);

% Fintech
x   = -0.75; %x = 0; % no lender frictions
par{11}=x;

q   = 0.05;  qmat = [0 q]';
pq  = 0.07; 
pprob1 =[1-pq pq];
qfmat = tensorfun(@plus,fmat,qmat); qfprob1 = kron(fprob1,pprob1); nqf = length(qfmat);
qfprob = repmat(qfprob1,nqf,1);


% Fintech + social effect
e  = 0.085; 
emat = [0 e]';
pe = 0.06;
eprob1 = [1-pe pe];
maxqe = [0 q e]'; maxqeprob = [(1-pq)*(1-pe) pq*(1-pe) pe];
qfemat = tensorfun(@plus,fmat,maxqe); qfeprob1 = kron(fprob1,maxqeprob); nqfe = length(qfemat);
qfeprob = repmat(qfeprob1,nqfe,1);


% initial state
T = 1;        % #. irf horizon
nhh = 100000;   % #. hh

distb = [20 20 40 14 3 2 1]/100;
BINI  = bmat(gendist(distb,1,nhh));    
ltv = 100*BINI/exp(median(pmat));

AINI  = 0.2*rand(nhh,1);AINI(1:0.4*nhh)=0;AINI(end-0.2*end:end)=0.25;


  % income
zini_M_ini  = genshock(zmat,ztran,muz,sz,nhh,1,0);  
zini_M = genshock(zmat,ztran,muz,sz,nhh,T,zini_M_ini);

distz_H = [4.81 10.8 20.28 21.16 24.04 13.08 5.83]/100;
zini_H_ini = zmat(gendist(distz_H,1,nhh));
zini_H = genshock(zmat,ztran,muz,sz,nhh,T,zini_H_ini);

  % fintech shocks
fini  = fmat(gendist(fprob(1,:),1,nhh));
qini  = qmat(gendist([1-pq pq],1,nhh));
qfini = fini+qini;

eini   = emat(gendist(eprob1,1,nhh));
qfeini = fini+max(qini,eini);
